Dynamics of shallow dark solitons in a trapped gas of impenetrable bosons 
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The dynamics of linear and nonlinear excitations in a Bose gas in the Tonks-Girardeau (TG) 
regime with longitudinal confinement are studied within a mean field theory of quintic nonlinearity. 
A reductive perturbation method is used to demonstrate that the dynamics of shallow dark soli- 
tons, in the presence of an external potential, can effectively be described by a variable-coefficient 
Korteweg-de Vries equation. The soliton oscillation frequency is analytically obtained to be equal 
to the axial trap frequency, in agreement with numerical predictions obtained by Busch et al. [J. 
Phys. B 36, 2553 (2003)] via the Bose-Fermi mapping. We obtain analytical expressions for the 
evolution of both soliton and emitted radiation (sound) profiles. 



Dark solitons are fundamental nonlinear excitations of 
the nonlinear Schrodinger (NLS) equation with repul- 
sive interactions, which have attracted much attention 
in the field of atomic Bose-Einstein condensates (BECs), 
where they have been studied experimentally [1] and the- 
oretically in the framework of the Gross-Pitaevskii (GP) 
equation [2-9]. In particular, under conditions of tight 
transverse confinement, for which the gas becomes quasi- 
one-dimensional, dark solitons in longitudinal harmonic 
traps have been shown to oscillate with frequency 0/v2 f 
where is the axial trap frequency [2-4]. Thermal [5] 
and dynamical [6] instabilities, including sound emission 
phenomena [7,8], have also been investigated. Further- 
more, apart from the traditional rectilinear solitons, ring- 
shaped dark solitons, have also been proposed [9]. 

For very tight transverse confinement and sufficiently 
low densities (n <C l/£, where £ the healing length of 
the system) [10], an atomic gas with repulsive interac- 
tions enters the Tonks-Girardeau (TG) regime [11], where 
it behaves like a one-dimensional gas of impenetrable 
bosons. This regime, which is currently under intense 
experimental investigation [12], has attracted consider- 
able theoretical attention both in the absence [13-15], 
as well as in the presence of axial confinement [16-18]. 
Dark solitons have also been predicted to occur in the TG 
regime, either upon employing the Bose-Fermi mapping 
theorem [14,19], or via a mean-field approach of quintic 
nonlinearity arising in this limit [20,21] (see also relevant 
work for fermionic systems [22]). However, some impor- 
tant issues of soliton dynamics, such as its instability to 
sound emission due to the longitudinal confinement, and 
the profiles of the emitted sound waves have not been 
discussed in detail. It is the aim of this paper to address 
these points analytically in the case of shallow solitons. 

Our analysis is based on the modified NLS equation 
with a quintic nonlinearity (as opposed to the cubic non- 
linear GP equation valid in 3D systems [23]) given by 
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Here ifi is the order parameter of the system (normalized 



to the number of particles N), m the atomic mass and lu x 
the axial confining frequency. The above long wavelength 
equation [20] , originally derived in [24] from a renormal- 
ization group approach, has also been obtained by other 
techniques, such as (i) the Kohn-Sham density functional 
theory for bosons [21,25,26], (ii) an in-depth study of 
the energy dependence of the many-body effective inter- 
action [27], (iii) a rigorous analysis of the many-body 
Schrodinger equation [28]. This equation is expected to 
be valid for weak density modulations, as long as the 
number of atoms exceeds a certain minimum value (typ- 
ically much larger than 10), for which oscillations in the 
density profiles become essentially suppressed [20,27,29]. 
Note that an alternative density-functional approach to 
deal with the excitations and dynamics in the TG limit 
has been recently proposed [30]. 

The validity of Eq. (1) has been criticized in certain 
regimes where the outcome is sensitive to the exact phase 
of the order parameter, e.g., interference experiments on 
a torus, where effects beyond the realm of mean field 
theory arise, and for which Eq. (1) overestimates the co- 
herence properties [31]. The latter treatment, however, 
ignores longitudinal confinement which is the key source 
of the dynamics studied in this paper. An alternative 
exact treatment [19] including the harmonic confinement 
has predicted soliton oscillations at the trap frequency 
(this result was explicitly demonstrated numerically for 
deep dark solitons). This paper complements the result 
of [19], in that the same oscillation frequency is obtained 
analytically for shallow dark solitons, for which the em- 
ployed hydrodynamic approximation is expected to be 
valid. 

Solitons moving on a background density gradient are 
prone to sound emission [2-4,8], but are stabilized against 
decay in a single harmonic trap by the continuous sound 
emission-reabsorption cycles [7] . We derive analytical re- 
sults for both the soliton amplitude and speed as well 
as associated sound profiles. Starting from Eq. (1), 
our analytical predictions are obtained by means of a 
reductive multi-scale perturbative method, which yields 
a variable coefficient Korteweg-de Vries (KdV) equation 
for the dynamics of shallow dark solitons. Note that the 



1 



above quintic NLS equation is a model of mathemati- 
cal interest in its own right, as "traditional" analytical 
techniques, such as the adiabatic perturbation theory for 
dark solitons in BECs [4,9], cannot be applied for such 
nonlincarities. However, a similar technique as the one 
presented here has been succesfully applied to quasi-one- 
dimensional condensates in the GP limit [3] . To start our 
analysis, we express Eq. (1) in dimensionless form, 

iu t = ~l;u xx + |u| 4 u + V(x)u, (2) 

where the subscripts denote partial differentiation, while 
the variables t, x and the function u are measured, re- 
spectively, in units of m/(7i7r 2 n 2 ), l/(nn ) and 2 1 / 4 y / n " 
(where n Q = \tp \ 2 the peak density of the gas). The ex- 
ternal potential in Eq. (2) is V(x) = (l/2)(Trn a x )~ 4: x 2 , 
where a x = (muj x ) is the harmonic oscillator length 
in the axial direction. Taking into regard that the pa- 
rameter (irn a x )~ is apparently small, it is convenient 
to define the small parameter e = il~ 2 ^ 3 (TTn a x )~ 4 ^ 3 [Q is 
a parameter of order O(l)], which will be used in the per- 
turbation analysis to be presented below. This way, the 
external potential takes the form V(X) = (l/2)f2 2 A 2 , 
i.e., it is a function of the slow variable X — e 3 ' 2 x, while 
f2 expresses the strength of the magnetic trap, or the 
normalized axial trap frequency. By analogy to [29], we 
introduce the Madelung transformation u = y/nexp(i(/>) 
(n and <f> denote the density and phase respectively) to 
reduce Eq. (2) to the following hydrodynamic equations, 

n t + {n4> x ) x = 0, (3) 

<k + \€ + n 2 - ^n-V V /2 )** + V(X) = 0, (4) 

which arc similar to the ones that have been employed to 
discuss the crossover from TG to BEC regime [17,18,32]. 
The ground state of the system may be obtained upon 
assuming that the atomic velocity v = (j> x = (i.e., no 
flow in the system) and <p t = — /i (dimensionless chem- 
ical potential). Then, as Eq. (3) implies that n = uq 
is time-independent in the ground state, we assume that 
no = no(X). As the quantum pressure term in Eq. (4) is 
of order 0(e 3 ), to leading order in e [to 0(1)], we obtain 

n (X) = V"o - V(X), (5) 

in the region where u > V(X), and n = outside. 
Equation (5) gives the density profile in the so-called 
Thomas- Fermi (TF) approximation [23]. It follows from 
Eq. (5) that in the case of the harmonic trap the ax- 
ial size of the gas is 2L, where L = y/2(i /Q- We now 
consider the propagation of small-amplitude linear exci- 
tations (e.g., sound waves) of the ground state, by seeking 
solutions of Eqs. (3)-(4) of the form, 

n = n (A) + eh(x, t), (j) = — (tot + e<f>(x, t), (6) 

where the functions n and <j> describe the linear excita- 
tions. Substituting Eqs. (6) in Eqs. (3)-(4), to order 



0(1) we recover the TF approximation, while to order 
0(e) we obtain the following system of linear equations, 

n t + n (f> xx =0, 4> t + Znon - -j—n X x = 0. (7) 

The dispersion relation of the excitations can easily be 
obtained, upon considering plane-wave solutions of Eqs. 
(7) of the form (n, cfi) — (ho, <j>o) exp[i(kx — wt)], where 
n and (j) are independent of x and t (but may depend 
on X). This way, we readily obtain the equation, 



ui = ± ] j2n 2 (X)k 2 + Jfc 4 (8) 

which is a Bogoliubov-type excitation spectrum, but with 
the excitation frequency u being a function of the slow 
variable X. The speed of sound is local and is given by 

c s = sV2n Q (X), (9) 

where s = sign(cg) = ±1, i.e., the sound may propagate 
in two opposite directions. Note that the local character 
of the speed of sound is due to the presence of the ex- 
ternal potential, which bears resemblance to the sound 
propagation in slowly- varying nonuniform media [33]. 

We now employ the reductive perturbation method 
[34,35] (see also [3] for a relevant study in BECs) to ex- 
amine the evolution of the nonlinear excitations (e.g., 
solitons), of the ground state. As Eqs. (3)-(4) are in- 
homogeneous, we introduce a new slow time-variable, 
T = e 1 / 2 (t - f C- X {rf)dxf), where C is the (local) ve- 
locity of the nonlinear excitations, to be determined in a 
self-consistent manner. Also, we introduce the following 
asymptotic expansions for the functions n and <f), 

n = n (A) + em (A, T) + e 2 n 2 (A, T) -\ , 

= -(L t + e 1/2 ^i (A, T) + e 3 / 2 2 (A, T) -\ . (10) 

Substituting the expansions (10) into Eqs. (3)- (4), we 
obtain the following results: First, to order 0(1), Eq. 
(4) leads to the TF approximation [see Eq. (5)]. Then, 
to orders O(e) and 0(e 3 / 2 ), Eqs. (4) and (3) respectively 
lead to the following system of linear equations: 

4>it + 2n ni = 0, niT + C~ 2 n (f>iTT = 0. (11) 

The compatibility condition of Eqs. (11) is the equation 
1 — 2noO~ 2 = 0, which determines the unknown velocity 

0, which, actually, is the same as the speed of sound, 

1. e., C = c s [c.f. Eq. (9)]. On the other hand, Eqs. (11) 
lead to the following equation, 

MX,T) = -2n (X) J ni (X,T')dT', (12) 

connecting the amplitude n\ and the phase <j>i. 

Proceeding to the next order, namely to order 0(e 2 ) 
and to order 0(e 5 / 2 ), Eqs. (4) and (3) respectively read: 
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/>2T + 2n n 2 



-C- 2 4>\ T + -C- 2 n^n 1TT , (13) 



n 2 T + C 2 n cj>2TT = -C 2 {nicj>i T ) T 



+C 1 (n 0X 4>iT + 2n 4> 1X T) + 
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-n (pi T . 



(14) 



The compatibility condition of Eqs. (13) and (14) yields 
the equation 1 — 2n 2 l C~ 2 = for the velocity C, along 
with the following nonlinear evolution equation for m, 



nix 



C 

—^n-invr 
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n lTTT 



dX 



(lnno)n 1; (15) 



which is obtained upon also using Eq. (12). Equation 
(15) has the form of a KdV equation with variable coef- 
ficients, which has been used to describe shallow water- 
waves over variable depth, or ion-acoustic solitons in in- 
homogeneous plasmas [35]. The inhomogeneity- induced 
dynamics and dissipation of the KdV solitons has been 
studied analytically [36,37]. Below, we will employ these 
results to analyze the coherent evolution and dissipation 
of dark solitons in an atomic gas in the TG limit. 

Introducing the transformations \ — J (32riQC)~ 1 dX , 
t = T and n\ = (3/32nl)v(x, t), we put Eq. (15) into 
the form, 



6vv T + v T 



(16) 



where A(x) = 2(lnno) x . In the case A = 0, i.e., for a 
homogeneous gas with n a (X) = n Q = const., Eq. (16) is 
the completely integrable KdV equation, which possesses 
a single-soliton solution of the following form [38] , 



-2ft 2 scclr Z, Z 



C(X)} 



(17) 



where C(x) = 4k 2 x + Co is the soliton center (with 
d(/dx = 4k 2 being the soliton velocity in the t-\ refer- 
ence frame), while n and Co are arbitrary constants pre- 
senting the soliton's amplitude (as well as inverse tempo- 
ral width) and initial position respectively. Apparently, 
Eq. (17) describes a density notch on the backround 
density n 0l with a phase jump across it [see Eq. (12), 
which implies that 4>i ~ tanhZ] and, thus, it represents 
an approximate dark soliton solution of Eq. (2). 

In the general case of the inhomogeneous gas (i.e., in 
the presence of the trapping potential), soliton dynam- 
ics can still be studied analytically, provided that the 
right-hand side of Eq. (16) can be treated as a small 
perturbation. As A(x) is directly proportional to the 
density gradient, such a perturbative study is relevant 
in regions of small density gradients, which is consistent 
with the use of the local density approximation. In this 
case, employing the perturbation theory for solitons [39], 
we express the solution of Eq. (16) as, 



v = v s +vr, 



(18) 



where vs is the soliton part, which has the same func- 
tional form as in the unperturbed homogeneous case (c.f. 



Eq. (17)), but with the soliton parameters k and £ being 
now unknown functions of %. The contribution vr, being 
of the same order of smallness as A, denotes the radia- 
tion part of the solution (i.e., the sound profile) due to 
the effect of axial inhomogeneity. Following [37], we first 
derive the following evolution equations for the soliton's 
amplitude k(x) and center C(x): 



dn _ 2 ^ d( _ Ak 2 + x 
dx 3 dx 3k 



(19) 



These equations can be solved analytically and the re- 
sults, expressed in terms of the slow variable X, read: 



(20) 



K 2 (0) 



8\/2srio /3 (0) ^o 



X 



n (X) 



n - 7/3 (x')dx' 



3k(0) 



MO) 



-1/3' 



(21) 



where n(0) and Uq(0) are respectively the soliton ampli- 
tude and density at X = 0, while s = ±1 represents the 
two possible directions of the soliton propagation. Addi- 
tionally, we find the following approximate (for \Z\ 3> 1) 
expression for the radiation part of the solution (i.e., 
sound wave emitted by the soliton), 



32V2 S ^ /3 (Q) 8/3 dnojX) 



VR W - 



3/t(0) 



(xy 



dX 



-(1-tanhZ), (22) 



where Z = k(X) [t — ((X)] is the soliton phase, with 
k(X) and ((X) given by Eqs. (20)-(21). 

Based on the above results for the evolution of the soli- 
ton parameters, we will now show that the dark soliton 
will display an oscillatory motion in the harmonic trap 
V{X) = (l/2)n 2 X 2 . This can be done upon deriving the 
phase of the soliton, which, to order 0(e 3 / 2 ), reads 



Z = £ 1/2 ^o /3 k(0) 



1 - 



r dx 
J ~c 



X(0)Mo 



4 (X 
3 

2/3 



(23) 



(recall that L = v / 27*o/^ defines the axial size of the 
gas). Then, looking along the characteristic lines of soli- 
ton motion, it is easy to show that the position of the 
soliton satisfies the following equation of motion, 



dX_ 
~dt 



8V2snLn (X) 



80L + ev / 2K 2 (0)/x 2/3 n a (X) ' 



(24) 



For sufficiently small e the second term in the denomina- 
tor can be neglected and the separable resulting equation 
can readily be integrated. In particular, taking into re- 
gard Eq. (5) for a parabolic trap, we find that 
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X = Lsin(ftt), 



(25) 



which demonstrates that a shallow dark soliton in the 
TG limit oscillates at the trap frequency. 

In conclusion, we have developed a systematic analyti- 
cal approach, based on a reductive perturbation method, 
to study the linear and nonlinear excitations of a Bose 
gas of impenetrable bosons (i.e., in the Tonks-Girardeau 
limit). We have recovered the Bogoliubov spectrum of 
linear excitations, with excitation frequencies (and speed 
of sound) varying slowly along the axial direction. Ad- 
ditionally, we have shown that shallow dark solitons can 
be described by an effective Korteweg-de Vries equation 
with variable coefficients. We have found analytical ex- 
pressions for the inhomogeneity-induced evolution of the 
soliton parameters (amplitude, width, position, velocity) 
and the profile of the sound emitted by the soliton. 

Our results are based on the quintic nonlinear 
Schrodinger equation [20,21,24-29], which is expected to 
be valid for weak density modulations. This approach 
enables analytical results to be obtained, and the oscil- 
lation frequency of shallow solitons is thus found to be 
identical to the one obtained via numerical simulations 
based on the Bose- Fermi mapping [19]. This result seems 
to additionally justify a posteriori the use of Eq. (1) for 
shallow solitons. In particular, our work demonstrates 
the validity of the reductive perturbation method for dy- 
namics of shallow solitons in the Tonks-Girardeau (TG) 
regime, while an earlier study [3] confirmed its useful- 
ness in the opposite regime of Bose-Einstein condensation 
(BEC). The crossover between these two regimes [17] is 
important, and diagnostics (e.g. [32]) for the degree of 
"impenetrability" of a trapped one-dimensional Bose gas 
are required for interpreting current and future experi- 
ments [12]. The presented technique paves the way for 
investigating the crossover in the oscillation frequency 
of shallow dark solitons [ranging from Q/V2 in the BEC 
regime, to Q in the TG regime (f2 is the trap frequency)], 
upon suitably generalizing the nonlinearity to a smooth 
function interpolating between these two regimes [17,30]. 
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